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Abstract —In this paper, in order to enhance the numerical 
stability of the unscented Kalman filter (UKF) used for power 
system dynamic state estimation, a new UKF with guaranteed 
positive semidifinite estimation error covariance (UKF-GPS) is 
proposed and compared with five existing approaches. Including 
UKF-schol, UKF-k, UKF-modlfied, UKF-AQ, and the square- 
root unscented Kalman filter (SR-UKF). These methods and the 
extended Kalman filter (EKF) are tested by performing dynamic 
state estimation on WSCC 3-machine 9-bus system and NPCC 48- 
machine 140-bus system. For WSCC system, all methods obtain 
good estimates. However, for NPCC system, both EKF and the 
classic UKF fall. It is found that UKF-schol, UKF-k, and UKF- 
AQ do not work well in some estimations while UKF-GPS works 
well in most cases. UKF-modifled and SR-UKF can always work 
well, indicating their better scalability mainly due to the enhanced 
numerical stability. 

Index Terms —Extended Kalman filter, dynamic state esti¬ 
mation, nonlinear filters, nonlocal sampling effect, numerical 
stability, phasor measurement unit (PMU), positive semidefinite, 
square-root unscented Kalman filter, synchrophasor, unscented 
Kalman filter. 
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Cross covariance of the state and measure¬ 
ment. 

Process noise and measurement noise col¬ 
umn vectors. 

Constant covariance matrices of q and r. 
Cholesky factor (matrix square root) of the 
estimation error covariance P. 

Weights for the mean and the covariance of 
the state or measurement. 

Column vector of the states. 

Sigma points and predicted sigma points. 
Column vector of the measurements. 
Predicted measurement. 

Propagated sigma points by the measure¬ 
ment function. 

Rotor angle in rad. 

Rotor speed and rated rotor speed in rad/s. 
Voltage source. 

Column vectors of all generators’ real and 
imaginary parts of the voltage source on 
system reference frame. 

Internal field voltage in pu. 

Terminal voltage phasor. 

Terminal voltage at q and d axes in pu. 
Transient voltage at q and d axes in pu. 
Real and imaginary part of the terminal 
voltage phasor. 

System state error averaged for one type of 
state (6, oj, e^, or e^j) over a time period. 
Number of generators. 

Number of generators with classical model 
and fourth-order transient model. 

Number of PMUs. 

Set of generators with second-order classi¬ 
cal model and fourth-order transient model. 
Set of generators where PMUs are installed. 
Generator inertia constant in second. 
Terminal current phasor. 

Current at q and d axes in pu. 

Real and imaginary part of the terminal 
current phasor in pu. 

Damping factor in pu. 

Number of states, inputs, and outputs. 
Electrical active output power in pu. 
System and generator base MVA. 
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Mechanical torque and electric air-gap 
torque in pu. 

Open-circuit time constants for q and d axes 
in second. 

Synchronous reactance at q and d axes in 
pu. 

Transient reactance at q and d axes in pu. 
Admittance matrix of the reduced network 
only consisting of generator^ 

The Ah row of Y. 

Cholesky factor of a matrix. 

Rank 1 update to Cholesky factorization. 
Obtain the eigenvalue and eigenvector of a 
matrix. 

Create diagonal matrix or get diagonal el¬ 
ements of matrix. 

Orthogonal-triangular decomposition of a 
matrix. 

Matrix square root of a positive semidefi- 
nite matrix P, which is a matrix S — VP 
such that P = SS^. 

Real part and imaginary part. 

The Ah column of a matrix. 

Columns of a matrix belonging to a set S. 
Frobenius norm of a matrix. 

2-norm of a vector. 

Elementwise product and matrix product. 


I. Introduction 

S TATE estimation is an important application of the energy 
management system (EMS). However, the widely studied 
static state estimation 0-0 assumes that the power system 
operates in quasi-steady state, based on which the static 
states of the system, i.e. the voltage magnitude and phase 
angles of all buses, are estimated by making use of the 
Supervisory Control and Data Acquisition (SCADA) and/or 
phasor measurement unit (PMU) measurements. 

Static state estimation is important for power system moni¬ 
toring and also provides input data for other important applica¬ 
tions in EMS, such as automatic generation control (AGC) and 
optimal power flow (OPE). However, it may not be sufficient 
for good system monitoring and situational awareness as the 
power system becomes more dynamic due to the increasing 
penetration of renewable generation that has very high un¬ 
certainty and variation. Therefore, accurate dynamic states of 
the system obtained from real-time dynamic state estimation 
(DSE) facilitated by high-level PMU deployment has thus 
become essential. With the high global positioning system 
(GPS) synchronization accuracy, PMUs can provide highly 
synchronized measurements of voltage and current phasors in 
high sampling rate, thus playing a critical role in achieving 
real-time wide-area monitoring, protection, and control. 

Power system DSE has been implemented by different 
types of Kalman filters. The most common application of the 
Kalman filter (KE) 0 to nonlinear systems is in the form 

*The elements of Y are constant if the difference between x'^ and x'^ is 
ignored ( JTj). 


of extended Kalman filter (EKE) fTOl , | [TT| , which linearizes 
all nonlinear transformations and substitutes Jacobian matrices 
for the linear transformations in KE equations, based on the 
assumption that all transformations are quasi-linear. Power 
system DSE has been implemented by EKF 1121, | [T3| . 

Although EKF maintains the elegant and computationally 
efficient recursive update form of the KE it works well only 
in a ‘mild’ nonlinear environment due to the first-order Taylor 
series approximation for nonlinear functions m- It is sub- 
optimal and can easily lead to divergence. The linearized 
transformations are reliable only when the error propagation 
can be well approximated by a linear function. Also, the 
linearization can be applied only if the Jacobian matrix exists. 
Even if the Jocobian matrix esists, calculating it can be a 
difficult and error-prone. 

The unscented transformation (UT) GD was developed to 
address the deficiencies of linearization by providing a more 
direct and explicit mechanism for transforming mean and 
covariance information. Based on UT, Julier et al. @-1111 
proposed the unscented Kalman filter (UKE) as a derivative- 
free alternative to EKE in the framework of state estimation. 
The UKE has been applied to power system DSE, for which 
no linearization or calculation of Jacobian matrices is needed 
1191, |20|. However, in | [T0 and pO) UKE is only applied to 
estimate the dynamic states for the single-machine infinite-bus 
system or WSCC 3-machine system. 

It is not surprising that UKE has not been applied to larger 
power systems. As has been pointed out in and pT[ , both 
EKE and UKE can suffer from the curse of dimensionality 
and the effect of dimensionality may become detrimental in 
high-dimensional state-space models with state-vectors of size 
twenty or more, especially when there are high degree of 
nonlinearities in the equations that describe the state-space 
model, which is exactly the case for power systems. 

Therefore, even if classic UKE has good performance for 
small systems, it might not work at all for large power systems. 
We will show that it is the numerical stability that mainly 
limits the scalability of the classic UKF Specifically, when 
the estimation error covariance is propagated, it sometimes 
cannot maintain the positive semidefiniteness, thus making its 
square-root unable to be calculated. 

In this paper, we introduce and compare six techniques that 
can be used to enhance the numerical stability of UKE, includ¬ 
ing the EKE/UKE toolbox approach | [22|, properly setting a 
parameter of unscented transformation |23| , the modified UKE 
approach GZ)^ adding an extra positive definite matrix p^ , 
p5| , the UKE with guaranteed positive semidifinite estimation 
error covariance (UKE-GPS) proposed in this paper, and the 
square-root UKE (SR-UKE) @. 

The remainder of this paper is organized as follows. Section 
[n| briefly introduces the unscented transformation and the 
classic UKE procedure. Sectionj^discusses six techniques for 
enhancing the numerical stability of the classic UKE. Section 
[^explains how Kalman filters can be implemented for power 
system dynamic state estimation. Section |V] tests the proposed 
methods on the WSCC 3-machine 9-bus system and NPCC 
48-machine 140-bus system. Finally the conclusion is drawn 
in Section FVll 
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II. Unscented Kalman Filter 
A discrete-time nonlinear system can be described as 

r Xk = f{xk-i,Uk-i) + Qk-I (la) 

\ yk = Hxk,Uk) -bT-fe, (lb) 

where Xk G M", ttfe S K", and S are, respectively, state 
variables, inputs, and observed measurements at time step fc; 
the estimated mean and estimation error covariance are m and 
P; f and h are vectors consisting of nonlinear state transition 
functions and measurement functions; qi^_i ^ N{0,Qi._i) is 
the Gaussian process noise at time step A: — 1; ~ N{0, Rk) 

is the Gaussian measurement noise at time step k\ and Qk_i 
and Rk are covariances of Qk-i and rk- 


A. Unscented Transformation 


Unscented Transformation (UT) is proposed based on the 
idea that “it is easier to approximate a probability distribution 
than it is to approximate an arbitrary nonlinear function or 
transformation” GD- A set of sigma points are chosen so that 
their mean and covariance are m and P. The nonlinear func¬ 
tion is applied to each point to yield a cloud of transformed 
points and the statistics of the transformed points can then be 
calculated to form an estimate of the nonlinearly transformed 
mean and covariance. 

Specifically, a total of 2 n -b 1 sigma points (denoted by X) 
are calculated from the columns of the matrix rjs/P as 


f = m 


1 X^''^ = m + 

rjs/P 

[ AfO) = m - 

' 

, 


with weights 


i = 1,..., n 
i = n -b 1,..., 2 n 


' ws 

irg 


A 


n -b A 
A 

n -b A 


+ (1 


1 

2(n -b A) ’ 
1 

2(n -b A) ’ 


-«'+/?) 

i = 1,..., 2 n 

i = 1,..., 2 n, 


(2a) 

(2b) 

(2c) 

(3a) 

(3b) 

(3c) 

(3d) 


where the matrix square root of a positive semidefinite matrix 
P is a matrix S = s/P such that P = SS^, Wm and Wc 
are respectively weights for the mean and the covariance, rj = 
s/n + X, A is a scaling parameter defined as A = a‘^{n+K)—n, 
and a, /3, and k are constants and a and (3 are nonnegative. 


B. Unscented Kalman Filter 

Assume the initial estimated mean and the initial estimation 
error covariance are mg and Pg, UKF can be performed in a 
prediction step and an update step, as in Algorithms and 


III. Unscented Kalman Filter with Enhanced 
Numerical Stability 


Here, we propose a UKF-GPS method (see Section III-Ei 
and introduce five other approaches to enhance the numerical 


Algorithm 1 UKE Algorithm; Prediction Step 
I: calculate sigma points 

Xk-i = [ ruk-i ■ ■ -mfc.ij 
2n+l 

+ V [Ori,l s/ P k-1 — sj Pfe-l] ■ (4) 

2 : evaluate the sigma points with the dynamic model func¬ 
tion 

Xk = f{Xk-i). (5) 


3: estimate the predicted state mean 




(6) 


2=0 


4: estimate the predicted error covariance 


2n 


Pk^Y. - mfKX^^k - m-)T + Qk_,. (7) 


i=0 


5: calculate the predicted sigma points 

0n,l 


Pk -\^Pk 


■ ( 8 ) 


2n+l 


6: evaluate the propagated sigma points with measurement 
function 

yk = h{X-k). (9) 

7: estimate the predicted measurement 




( 10 ) 


i=0 


Algorithm 2 UKF Algorithm: Update Step 
I: estimate the innovation covariance matrix 

2n 

Py.y. = E fc - yfc ) fc -ykY +Ru. (11) 

2 = 0 

2: estimate the cross-covariance matrix 

2n _ 

P^.y. = E i^lk - mf) {y-k -y'k,^ 


2 = 0 

3: calculate the Kalman gain 

_ p P~1 

4: estimate the updated state 

ruk = mf + Kk{yk-yk)- 
5: estimate the updated error covariance 

Pk = Pk ~ R-kPy^y^Kk ■ 


( 12 ) 

(13) 

(14) 

(15) 


stability of the classic UKF. We also summarize and discuss 
the advantages and disadvantages of these approaches. 

In Section [ri-B[ the estimation error covariance in Algorithm 
should be positive semidefinite, because its square root is 
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required in order to obtain the sigma points, as shown in Q 
and (§. However, through propagation the estimation error 
covariance can lose positive semidehniteness. 

As for why the estimation error covariance can lose positive 
semidehniteness for the classic UKF, it has been shown in 
the Appendix III of llz) that when k, a parameter used for 
unscented transformation, is negative it is possible to calculate 
a nonpositive semidehnite estimation error covariance. As 
mentioned in p7| , this problem is not uncommon for methods 
that approximate higher order moments or probability density 
distributions, as those described in fTO) , and | |28l . 

In flTt a useful heuristic is proposed as n + k = 3 which 
can minimize the moments of the standard Gaussian and the 
sigma points up to the fourth order. From Q it is seen that 
the distance of the sigma point from the mean is proportional 
to 77 = If the UKF procedure follows the heuristic 

n + K = 3, the desired dimensional invariance is achieved by 
canceling the effect of the system dimension n, thus avoiding 
the sampling of nonlocal effects that can lead to signihcant 
difficulties in worst cases p3) , p^ . However, for a high 
dimension system with big n, the weight of the center point 

= ^ + {l-a^ + l3) 
n + X 

= 2-0^+ (16) 

can be negative. For a typical selection a = 1, /3 = 0, = 

\—n/3. When n > 3, will be negative and the calculated 
covariance may become nonpositive semidehnite. 


A. EKF/UKF Toolbox Approach 

In EKF/UKF toolbox p^ , when Pk-i or is not 
positive semidehnite, the function ‘schol’, which calculates 
the lower triangular Cholesky factor of a matrix, can still give 
an output. The ‘schol’ algorithm can be summarized as 


k=l 


Sn = 


0 , 


if s > e 
otherwise 


j-i 

ij — ^ ^ ^ik^jk 

k=l 


Srj = 


if Sjj > e 
otherwise, 


(17) 

(18) 

(19) 

( 20 ) 


where P is the covariance matrix, S is the output of the ‘schol’ 
function, and e = 2.22 x 10 “^® is the distance from 1.0 to the 
next largest double-precision number in MATLAB. If a matrix 
P is positive semidehnite, ‘schol’ can obtain a S matrix such 
that P = SS^. When P is positive semidehnite, the ‘schol’ 
can still get a matrix S but P = SS^ cannot be satished. 
However, by using this S the sigma points can be calculated 
and the estimation by UKF can at least continue to proceed. 
This approach for enhancing the numerical stability is called 
“UKF-schor. 


B. Selection of k 

When K is negative it is possible to calculate a nonpositive 
semidehnite estimation error covariance GD- Therefore, in 
1 23 1 it is suggested to choose /c > 0 to guarantee the 


positive semidehniteness of the the covariance matrix. Since 
the specihc value of k is not critical, a good default choice is 
K = 0 J^ . This approach is named as “UKF-k”. 

When K = 0 , the distance of the sigma point from the 
mean is proportional to y/n. As n increases, the radius of 
the sphere that bounds all the sigma points also increases 
| |2^ . Even though the mean and covariance of the prior 
distribution are still captured correctly, it does so at the cost 
of possibly sampling nonlocal effects, which can lead to 
signihcant difficulties if the nonlinearities in question are very 
severe. Therefore, although selecting k — 0 addresses the 
numerical instability problem in UKE, it picks up the nonlocal 
sampling problem. 


C. Modified UKF 

In GD a useful heuristic is proposed as n + k = 3 which 
can minimize the moments of the standard Gaussian and the 
sigma points up to the fourth order. This means that for a 
system with n > 3, k will be negative. In order to avoid 
obtaining a nonpositive, semidehnite covariance when using a 
negative k, a modihed UKE is proposed in GZ) for which the 
predicted error covariance in Q and the innovation covariance 
matrix in ( 111 are evaluated about the projected mean as 


Pk=Y. - Xo,kV + Qk-i 


2 = 0 
2n 


= ^ - Xo,k){X,,k - + Qk-1 ( 21 ) 


2 = 1 


2n 

Py.y. = E {yik - y^k) {yik - y^^kV + Rt. 

2=0 

2n 

= E iyT,k - yo,k) ( t . a - yo,kV +Rk ( 22 ) 

2=1 

It is shown in GZl that the modihed form ensures positive 
semidehniteness, and, in the limit (n -f k) —> 0 , the modihed 
UKE is the same as that of the modihed, truncated second- 
order hlter p 8 ) . This approach is called “UKF-modified”. 

D. Adding AQ 

In p4) and | |25) , an extra positive dehnite matrix AQk is 
added to the predicted covariance matrix in (0 as a slight 
modihcation of the UKF to improve the stability of UKF. It 
is shown that the estimation error of the UKF is bounded if 
AQk is set properly and the stability of UKF is improved. 
However, the precision of the estimation can be decreased. 
This approach is called “UKF-AQ”. Specihcally, the predicted 
error covariance in 0 becomes 

2n 

2=0 
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where Qi~_i — Qk-i + In |24| a nonlinear system 

with linear measurement functions are considered and no 


method is provided to design the additional AQ^. while in [251 
a nonlinear system with nonlinear measurement functions are 
considered and a heuristic method is provided to design AQk- 


E. UKF-GPS 

If Pk-i or is nonpositive semidefinite, the UKF-GPS 
will execute the nearest symmetric positive definite (nearPD) 
algorithm (a R function in ‘Matrix’ package p0|), as shown 
in Algorithm by which a symmetric positive semidefinite 
matrix nearest to Pk-i or P^ in Frobenius norm can be ob¬ 
tained. The input Xq can be Pk-i or P^ and is converted to 
the output X, which guarantees the positive semidefiniteness 
and substitutes Pk-i or P^. 

The ‘nearPD’ algorithm adapts the modified alternating 
projections method in pT| and then adds procedures to force 
positive definiteness by ‘posdefify’ (a R function in ‘sfsmisc’ 
package) p^ , and to guarantee symmetric. The modified 
alternating projections method iteratively projects a matrix 
onto the set 5 = {K = e ; F > 0} by a 

modified interation due to Dykstra {AS is Dykstra’s 
correction), which incorporates a judiciously chosen correction 
to each projection that can be interpreted as a normal vector 
to the corresponding convex set pT) . As is mentioned in 
ID, general results in p4) and ]35| show that both X and 
Y converge to the desired nearest covariance matrix as the 
number of iterations approach infinity. The rate of convergence 
of Dykstra’s algorithm is linear when the sets are subspaces 
and the constant depends on the angle between the subspaces 
p^ . To force positive definiteness, the eigenvalues less than 
Eps are replaced by a positive value Eps. 

In Algorithm]^ ‘eig’ (eigen decomposition), ‘max’, ‘sqrt’ 
(square root), ‘diag’, ‘•’ (element-wise product), ‘x’ (matrix 
product), and ‘./’ (element-wise division) are MATLAB func¬ 
tions; V is the matrix of eigenvectors, d is the vector of 
eigenvalues; p is the elements that satisfy d > reigmax(d); 
\y]p is the columns of V that belong to p\ dp is the rows 
of d that belong to p; and ||A|| is the Frobenius norm, the 
matrix norm of an m x n matrix A with entry defined as 


lAII = 


\ 


EEi 

i=i i=i 


(24) 


E SR-UKF 

The calculation of the new set of sigma points at the predic¬ 
tion step requires taking a matrix square-root of the covariance 
matrix P by SS^ = P. For UKF, while the square-root of P 
is an integral part, it is actually still the full covariance P that 
is recursively updated. During the propagation, it is possible 
that P can lose its positive semidefiniteness. By contrast, 
in the implementation of SR-UKF, S is directly propagated, 
thus avoiding refactorizing P at each step. SR-UKF has been 
applied to power system DSE in p7)-p^. 

SR-UKF can be implemented by Algorithms and The 
filter is initialized by calculating the matrix square-root of the 


Algorithm 3 nearPD Algorithm 

1: 

initialize 



Let AS = On,n- 


2: 

modified alternating projections 



do 



Y = X 

(25) 


R=Y -AS 

(26) 


[V,d]^eigiR) 

(27) 


p ^ d> Teig max(d) 

(28) 


X = [V]p-[dp_^] X [V]l 

(29) 


n 

AS = X-R 

(30) 


while ||F-X||/||X|| >Teo„v 


3: 

guarantee positive definite 



[y,d]^eig(X) 

(31) 


Eps 3— Tposd max((i) 

(32) 


d{d < Eps) 3— Eps 

(33) 


diagX ^ diag(X) 

(34) 


X = Vdiag{d)V^ 

(35) 


D = sqrt [inax{Eps, diagX)./diag{X)) 

(36) 


X = diag(£>) X X-[D_^]. 

(37) 

4: 

guarantee symmetric 



0 

(38) 


estimation error covariance once via a Cholesky factorization 
as So = schol {Pq) where ‘schol’ is a function in EKF/UKF 
Toolbox that calculates the Cholesky factor of a matrix. 
The propagated and updated Cholesky factor is then used in 
subsequent iterations to directly form the sigma points. 


Correspondingly, (j5^-(|5^ in step 4 of Algorithmic replace 
the estimation error covariance update Q in Algorithm 
(57 1 -( 601 in Step 1 of Algorithm replace the innovation 
covariance update (111 in Algorithm |2 (|62| replaces 
for calculating Kalman gain; and (|64|-( 65]ireplace ([^ by 
applying p sequential Cholesky downdates to Sf^ where p is 
the number of outputs. 


In Algorithms [C and |C the ‘qr’ (orthogonal-triangular 
decomposition) and ‘cholupdate’ (Rank 1 update to Cholesky 
factorization) are MATLAB functions; ‘s’ denotes the sign of 
and will be ‘-I-’ if > 0 and ‘-’ otherwise. 


We first show why ([50|)-((53]l is equivalent to (0- Eor the 

' ' AG K3nx« 

]) , (39) 


matrix in (50i which is now denoted by A C as 

_n T 

A = ' ' ‘ 




2n 
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a QR decomposition can be performed as 

Ri 


A — QR — [0i Q 2 


O 2 


= Qi^i 


(40) 


are all 

T : 


where Q G G and <53 ^ 

unitary matrices (for a unitary matrix B, there is B' B — 
BB^ = I), Ri G I^nxra jg upper triangular matrix, 
QiRi is called the thin QR factorization or reduced QR 
factorization E3’ and there is 

Ri = Ri Qi QiRi — A~^ A 


= Y, - m~y + Q,_,. 


i—1 

(41) 

The Sf. in (51 1 on the left hand side of the arrow is actually 
Ri. Then for the on the left hand side of (52 1 we have 


HO) 


'RjRi + IWi°^l(Ao,k-m^), if > 0 

rJRi - - "ifc ), otherwise. 


(42) 


From ([4l]i-([42|, it is easy to obtain 

isiYsi 

2n 

=Pk 

By ( [ 5 ^ we convert the upper triangular matrix to a lower 
triangular matrix and for S'^ on the left side of (53 1 there is 


Qk-l 


(43) 


=P^. 


(44) 


As for why ([^-([hgi can replace ( [TT| l, it is similar to why 
([^-([^ is equivalent to (0 and thus will not be discussed in 
detail. The relationship between the obtained from (57 1 - 
( [60j l and the Py^y^ in ( [TT] i can be written as 

Sy.Sl^Py^y^ (45) 

and therefore the Kalman gain calculated by ( |6^ is equivalent 
to the one in ([T3]l. Then from (|6^-(|6^ we have 


(46) 


SjSk = {S^VS^ - iK>,Sy){K>,SyJ 


~ Pk ^kPyi^y^Kk 


Vk! 

T 


which is implemented by applying p sequential Cholesky 
downdates to where p is the number of outputs. Each 
Cholesky downdates uses one column of U as the column 
vector. Thus (64i-(65i is equivalent to (15 1 . 


G. Summary and Discussion 

The above-mentioned methods are summarized as follows. 

1) The UKF-schol approach does not solve the problem 
of the non-positive semi definiteness of the estimation 
error covariance but is only able to obtain an inaccurate 
Cholesky factor when the estimation error covariance is 
not positive semidefinite. 


Algorithm 4 SR-UKF Algorithm; Prediction Step 
1: calculate sigma points 

Xk-i — [ ruk-i ■ ■ ■ m,k-i ] 

2n+l 

-l-77[0„_i Sk-i -5fe-i]. (47) 

2 : evaluate sigma points with the dynamic model function 
Affc = /(Affc.i). (48) 


3: estimate the predicted state mean 

2n 

(49) 

i^O 

4: estimate the predicted square root of error covariance 





^ cholupdate ( 5 ^, (ATo.fe-m^), ‘s’) (52) 

SI ^ {SlY. (53) 


5: calculate predicted sigma points 


xi: = [ 


rriu • m, 


2n+l 


■]+r7[0„,i SI -SI]. (54) 


6: evaluate the propagated sigma points with measurement 
function 

yi = H^k)- (55) 


7: estimate the predicted measurement 

2n 


(56) 


i=0 


2) The UKF-k approach guarantees the positive semidefi¬ 
niteness of the estimation error covariance but discards 
the useful heuristic n -I- k = 3 for n > 3 and also picks 
up the nonlocal sampling problem. 

3) UKF-modified can also guarantee the positive semidef¬ 
initeness of the estimation error covariance. It is shown 
that under some conditions it is the same as that of the 
modified, truncated second-order filter in). 

4) For UKF-AQ approach, it is hard to select a proper 
extra positive definite matrix. The heuristic proposed 
in | [25) does not work for the case with non-positive 
semidefinite estimation error covariance. Also, if the 
process noise covariance is enlarged too much, the 
precision may be decreased; if it is not sufficiently 
enlarged, the estimation error covariance can still be 
non-positive semidefinite. It is more reasonable to find 
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Algorithm 5 SR-UKF Algorithm; Update Step 
1 : estimate the innovation covariance matrix 

[Q, ^ 

) (57) 


2n 






'Vk 


(58) 


cholupdate {yo.k-yk),‘^’) (59) 

(60) 


Svk^iSyj'. 

2: estimate the cross-covariance matrix 


P.kyk=J2^c\X-,-m-){yr^-y~f. (61) 


i=0 


(62) 


3: calculate the Kalman gain 

4: estimate the updated state 

mk=m]: + KkiVk-Vk)- ( 63 ) 

5: estimate the updated square root of error covariance 

u = KkSy^ (64) 

Sk = cholupdate(5'^, 17, (65) 


the nearest positive semidefinite matrix, as in UKF-GPS. 

5) UKF-GPS converts the estimation error covariance to 
the nearest positive semidefinite matrix whenever it loses 
positive semidefinateness. However, in some cases in or¬ 
der to guarantee positive semidehniteness the converted 
positive semidehnite matrix can be not so close to the 
original one, and may lead to decrease of precision. 

6 ) SR-UKF intrinsically guarantees the positive semidef¬ 
initeness of the estimation error covariance since the 
square root of the covariance rather than the covariance 
itself propagates. 

7) As for the implementation based on the classic UKF, 
UKF -K and UKF-AQ are easier than the others. UKF- 
schol needs to modify the Cholesky factor algorithm, 
UKF-modified needs to modify the covariance calcu¬ 
lation, and UKF-GPS requires to add the ‘nearPD’ 
algorithm. For SR-UKF, it does require more extensive 
changes of the Kalman hlter procedure. 

8 ) As for calculation efficiency, SR-UKF can be more 
efficient than other UKF-based methods, mainly because 
it makes use of powerful linear algebra techniques 
including the orthogonal-triangular decomposition and 
Cholesky factor updating. 

IV. Power System Dynamic State Estimation 

Here, we discuss how different Kalman filters are applied 
to dynamic state estimation. We apply the generator and 


measurement model in Section III.C of | |J7| , which can be 
used for multi-machine systems and allows both fourth-order 
transient generator model and second-order classical generator 
model. The terminal voltage phasor and terminal current pha- 
sor obtained from PMUs are used as the output measurements. 

Let ^4 and Q 2 respectively denote the set of generators 
with fourth-order and second-order model. The numbers of 
generators with fourth-order or second-order model, which 
are also the cardinality of the sets f /4 and Q 2 , are 54 and 
g 2 , respectively. Thus the number of states n = 4 p 4 -f 2 p 2 - 
For generator i G G4, the fast sub-transient dynamics and 
saturation effects are ignored and the generator model is 
described by the fourth-order differential equations in local 
d-q reference frame: 

Si =u}i- ujQ ( 66 a) 

Wi = (Tmi - Tsi - - Wo)) (66b) 

Zlii \ UJq / 

Sqi = 7 ^ {piii “ “ x'^i) idi) ( 66 c) 

dO* ^ ^ 

= + {x^i - x'^i) , ( 66 d) 

^qo* 


where i is the generator serial number. 

For generator i G Q 2 , the generator model is only described 
by the hrst two equations of (661 and the eC and are kept 
unchanged. The set of generators where PMUs are installed 
is denoted by Qp. For generator i G Gp, En = epi + jeu and 
Iti = *Ri + jhi can be measured and are used as outputs. T^i 
and Efii are used as inputs. 

The dynamic model (661 can be rewritten in a general state 
space form as 


r x = fAx,u) (67a) 

\ y = hc{x,u), (67b) 


where the state vector x, input vector u, and output vector y 
are respectively 


X = 

5^ 


T 

( 68 a) 

u = 




( 68 b) 

y = 

.4 

7 T VT] 

Ci Zr Zj j . 


( 68 c) 


The iqi, in, and Tn in (661 are actually functions of x\ 


T'Ri = e’n sin Si + cos Si (69a) 

T'li = e7 sin Si - e'n cos Si (69b) 

7, (69c) 

ipi = Re(/u) (69d) 

in = Im(/ti) (69e) 

iqi = (iu sin Si + ipi cos S^) (69f) 

in = («Ri sin Si - in cos Si) (69g) 

Cqi = Cqi - x'^iin (69h) 

Cdi = e'n + aiqjiq, (69i) 

Tel — Cqiiqi (69j) 
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Te, = (69k) 

In ( |69l l, the outputs ir and ii are written as functions of 
X. Similarly, the outputs e^i and eu can also be written as 
function of x: 


CR2 — sin Si Bqi COS Si 
Cli — Cqi sin Si Cdi COS Si. 


(70a) 

(70b) 


Note that we do not consider the dynamics of and 
Ef^ but assume they are constant and known, since the main 
objective of this paper is to discuss techniques that enhance 
the numerical stability of UKF. The dynamic state estimation 
with unknown inputs (T^ or Ef^) has already been discussed 
in and similar discussion under the framework of 

this paper will be specially investigated elsewhere. 

Similar to iz) and p2) , the continuous models in ( [66] l can 
be discretized into their discrete form as 

Xk = f{xk-i,Uk-i) (71a) 

yk = h{xk,Uk), (71b) 

where k denotes the time at kAt and the state transition 
functions / can be obtained by the modihed Euler method 
||43) as 


— ^k — 1 E f ci.^k — lj '^k — l)At 

'r ^ fc{^k,Uk) + f^{Xk-l,Uk-l) 
^ 2 
Xk = Xk-1 + fAt. 


(72) 

(73) 

(74) 


The model in (71 1 can be used to perform power system 
dynamic state estimation with different Kalman filters. 


V. Simulation Results 

Here, the UKF-GPS and SR-UKF are tested on WSCC 3- 
machine 9-bus system and NPCC 48-machine 140-bus system, 
which are extracted from Power System Toolbox (PST) |j44j. 
The EKF and classic UKF comes from EKF/UKF toolbox 
and the UKF-GPS and SR-UKF algorithms are implemented 
based on EKF/UKF toolbox. All tests are carried out on a 
3.2-GHz Intel(R) Core(TM) i7-4790S based desktop. 


A. Settings 

The simulation data is generated as follows. 

1) The simulation data is generated by the model presented 
in Section IV and the sampling rate is set to be 120 
samples per second. 

2) In order to generate dynamic response, a three-phase 
fault is applied at one bus of the branches with the 
highest line flows and is cleared at the near and remote 
end after 0.05 and 0.1 second. We do not consider the 
fault on lines either bus of which is a generator terminal 
bus because this can lead to the tripping of a generator. 

3) For each measurement, Gaussian noise with variance 
0.01^ is added. 

4) The sampling rate of the measurements is set to be 60 
frames per second to mimic the PMU sampling rate. 

5) Gaussian process noise is added and the corresponding 
process noice covariance is set as a diagonal matrix. 


whose diagonal entries are the square of 10% of the 
largest state changes, as in p^ . 

6) For WSCC system, one PMU is installed at the terminal 
bus of generator 3, and for NPCC system, 24 PMUs are 
installed at the terminal bus of generators 1, 2, 3, 4, 6, 
9, 10, 12, 13, 14, 16, 18, 19, 20, 21, 27, 28, 31, 32, 35, 
36, 38, 44, and 45; the PMU placements are determined 
by the method in which is based on maximizing 
the determinant of the empirical observability gramian. 

The considered Alters are set as follows. 

1) Dynamic state estimation is performed on the post¬ 
contingency system on time period [0,10 s], which starts 
from the fault clearing. 

2) The initial estimated mean of the system state is set to 
be the pre-contingency state. 

3) For all methods, a = 1 and /? = 0. For UKF-k method 
K = 0 and for all the other methods k = 3 — n. 

4) The initial estimation error covariance Pq is set as 


„2 r 

^ 9.9 

^9,94 

O9.94 

^ 9,9 

' ui-'-g 

^9,94 

O9.94 

O94.9 

^94,9 

^e ;-'94 

^ 94,94 

. O94.9 

^94,9 

f* 94.94 

1 _ 


(75) 


where rs and are chosen as O.Stt/ISO and 10 
and Te' and rp’ are set to be 10“^. 

As mentioned before, the covariance for the process 
noise is set as a diagonal matrix, whose diagonal entries 
are the square of 10% of the largest state changes | [42) . 
The covariance for the measurement noise is a diagonal 
matrix, whose diagonal entries are 0.01^, as in |42|. 
For UKF-AQ method, the additional positive definite 
matrix AQ is set to be 0.005^7„, as suggested in 
For ‘nearPD’, Teig = 10“® and Tconv = Tposd = 10” 

To quantitatively compare the estimation results, we define 
the following system state estimation error index 


5) 

6 ) 

7) 

8 ) 




g T, 

z 

i=l t=l 


E E 


gTs 


(76) 


where a; is a type of states and can be 5, oj, e^, or e^; xfl 
is the estimated state and is the corresponding true value 
for generator i at time step /; Tg is the number of time steps. 


B. WSCC 3-Machine System 

Different methods discussed in Section [111] are tested on the 
WSCC 3-machine system, as shown in Fig. [T| All genera¬ 
tors are assumed to have second-order classical model. The 
estimated state trajectories from different Kalman Alters are 
shown in Fig. for which a three-phase fault is applied at 
bus 8 of line 8 — 9, the line with the highest line flow. For 
this small system with only six states, there is no obvious 
numerical stability problem and all methods work well, even 
though for the UKF methods except the UKF-k method there 
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Fig. I. WSCC 3-machine 9-bus system. 
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UKF-modified Estimate 
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- -UKF-GPS Estimate 

- SR-UKF Estimate 
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375 ^ 



Time (second) 


Time (second) 


Fig. 2. Estimated states for WSCC 3-machine system. 


is K = — 3. In this case the estimation error covariance of UKF 
can keep its positive semidefiniteness during propagation. 

There are six branches no bus of which is a generator 
terminal bus. Since the three-phase fault can be applied to any 
one of the two buses, there are totally twelve possible fault 
scenarios. We perform DSE for each of them and calculate 
the average values of the system state estimation error index, 
which are listed in Table |I] The standard deviations of 
are also listed in the parentheses under ex- It is seen that 
all methods have small average error and standard deviation 
and among them SR-UKF has the smallest error and standard 
deviation. 

C. NPCC 48-Machine System 

As shown in Fig. the NPCC system | |44) represents the 
northeast region of the FI system. Twenty seven generators 
have fourth-order model and the other twenty one have second- 
order classical model. Thus there are a total of 150 states. 

We perform DSF for 50 times and for each of them a three- 
phase fault is applied at the from bus of one of the 50 branches 
with highest line flows. For all of the estimations, FKF fails to 
converge and the classic UKF encounters numerical stability 
problem because the estimation error covariance Pk-i or 
loses positive semidefiniteness at some time steps. Theoreti¬ 
cally, in this case the square root of Pk-i or P^ cannot 


TABLE I 

Average Estimation Error for WSCC 3-Machine System 


Filter 

es 


EKF 

0.0371 

(0.0167) 

0.394 

(0.0972) 

UKF-schol 

0.0526 

(0.0196) 

0.463 

(0.159) 

UKF-re 

0.0267 

(0.0141) 

0.306 

(0.102) 

UKF-modified 

omii 

(0.0148) 

0.312 

(0.104) 

UKF-AQ 

0.0464 

(0.0163) 

0.478 

(0.148) 

UKF-GPS 

0.0526 

(0.0196) 

0.463 

(0.159) 

SR-UKF 

0.0250 

(0.0136) 

0.295 

(0.0988) 



Eig. 3. Map of the NPCC 48-machine 140-bus system. The stai's indicates 
generators with classical model. 


be calculated. Thus the sigma points in Q or cannot be 
obtained and the estimation procedure has to halt. Note that 
both FKF and the classic UKF methods fail due to the infeasi¬ 
bility of the methods themselves rather than other factors such 
as the settings of the FKF/UKF toolbox or the convergence 
tolerance: 1) For both methods the FKF/UKF toolbox chooses 
typical parameters, and using these parameters both methods 
work well for the smaller WSCC 3-machine system but fail for 
the bigger NPCC 48-machine system for their poor scalability, 
which for FKF is because of the loss of nonlinear dynamics in 
the linearization of the nonlinear transformations and for the 
classic UKF is due to the above-mentioned numerical stability 
problem; 2) The estimated states from FKF quickly diverge 
to values with very large absolute values while the classic 
UKF cannot continue to perform estimation because of the 
numerical instability, and thus both methods fail not because 
of the choice of the convergence tolerance. 

The reason why the estimation error covariance can lose 
positive semidefiniteness for the classic UKF has been dis¬ 
cussed in Section Here we would like to emphasize that 
the selection of outputs or the measured values cannot cause 
the loss of positive semifefiniteness, since we use the same 
outputs and the same settings for simulation data generation 
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Time (Second) 


Fig. 5. Norm of relative error of the states for 11th estimation. 

and Kalman filters for both WSCC 3-machine system and 
NPCC 48-machine system and the estimation for WSCC 
system works very well. Also, the measurement configuration 
cannot be the cause since the numerical stability problem still 
exists even when all of the generators are installed with PMUs. 

In Fig. 1^ we show the estimation error index for each of 
the fifty estimations. We can see that UKF-schol, UKF-k, and 
UKF-AQ do not work well and can have very big estimation 
errors for several estimations. UKF-AQ even diverge for some 
estimations, for which the estimation error index is too big 
and thus is not shown. UKF-GPS works well for almost all 
estimations, except for the 10th estimation in which case it has 
smaller error for and than UKF-schol but has similarly 
big error of 6 and w. By contrast, UKF-modified and SR-UKF 
both work very well for all estimations due to their enhanced 
numerical stability and scalability. 

For the estimation error index of the rotor angle, the UKF- 
schol, UKF-k, UKF-AQ, and UKF-GPS get their maximum 
index among 50 estimations on the 11th, 36th, 3rd, 10th, 
respectively. In Figs. [5]-[^ we show the 2-norm of the relative 
estimation error of the states ||(a;f. — mf.)/a;fe ||2 where is 
the real states and mk is the estimated states. From these 
hgures it is seen that the UKF-schol, UKF-k, UKF-AQ, or 
UKF-GPS can get poor estimation while the UKF-modihed 
and SR-UKF can always guarantee much better estimation 
results. 

Similar to the WSCC system case, the average values of 
the estimation error index are also calculated, which are listed 
in Table It is seen that the average estimation error index 
and its standard deviation for UKF-modihed and SR-UKF are 
signihcantly smaller than the other methods. 

In the above estimations, we only apply three-phase faults to 
generate dynamic responses. To further validate the proposed 
approach, we now consider different types of faults, including 
three-phase fault, line to ground fault, line-to-line to ground 
fault, line-to-line fault, and loss of line. We perform DSE for 
50 times and for each of them a randomly selected type of 
fault is applied at the from bus of one of the 50 branches 
with highest line Hows. Similar to the case that only considers 



Fig. 6. Norm of relative error of the states for 36th estimation. 



Fig. 7. Norm of relative en'or of the states for 3rd estimation. 



Fig. 8. Norm of relative en'or of the states for 10th estimation. 
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Fig. 4. Estimation error index of the states by different methods for NPCC 48-machine system. 


TABLE II 

Average Estimation Error eor NPCC 48-Machine System 


Filter 

es 

Guj 

eg/ 

^e' 

^d 

EKF 

- 

- 


- 


34.018 

6.967 

0.0121 

0.027 

UKF-schol 


(191.176) 

(35.591) 

(0.035) 

(0.063) 


0.295 

0.406 

0.003 

0.015 

UKt-k: 


(0.498) 

(0.194) 

(0.002) 

(0.014) 


0.0145 

0.232 

0.002 

0.008 

UKF-modihed 


(0.004) 

(0.111) 

(0.0008) 

(0.007) 

UKF-AQ 

- 

- 

- 

- 


0.749 

0.519 

0.002 

0.010 

UKF-GPS 


(5.054) 

(1.027) 

(0.002) 

(0.008) 


0.017 

0.243 

0.002 

0.010 

SR-UKF 


(0.007) 

(0.132) 

(0.001) 

(0.010) 


TABLE III 

Average Estimation Error for NPCC 48-Machine System under 
Different Faults 


Filter 

65 


e-,. 

^e' 

®d 

EKF 

- 

- 


- 


4.268 

1.274 

0.007 

0.020 

UKF-schol 


(18.236) 

(3.745) 

(0.024) 

(0.060) 


0.174 

0.303 

0.002 

o.oio 

UKF-re 


(0.371) 

(0.212) 

(0.002) 

(0.011) 


0.013 

0.196 

0.001 

0.008 

UKF-modined 


(0.008) 

(0.170) 

(0.001) 

(0.010) 

UKF-AQ 

- 

- 

- 

- 


0.586 

0.546 

0.002 

0.010 

UKF-GPS 


(3.423) 

(0.903) 

(0.002) 

(0.010) 


0.013 

0.172 

0.001 

0.006 

SR-UKF 


(0.008) 

(0.139) 

(0.001) 

(0.008) 


three-phase faults, for all estimations EKF fails to converge 
and classic UKF encounters numerical stability problem. The 
average values of the estimation error index are listed in Table 
which shows that the UKF-modified and SR-UKF meth¬ 
ods have much better performance than the other methods. 
Compared with the case only considering three-phase faults, 
the estimation error is smaller, possibly because three-phase 
fault is the most severe fault and the corresponding dynamics 
can be farther away from normal operating conditions. 


As pointed out in |[T4) and EKF, UKF, and SR-UKF 
all have computational complexity of 0{n^). The average 
times for performing DSE by different Kalman filters are 
listed in Table HYi Here we list the calculation times for 
both only considering three-phase fault and randomly choosing 
different types of faults. Note that the time reported here is 
from MATFAB implementations and is not fully optimized. 
It can be greatly reduced by more efficient, such as C- 


based, implementations and by further optimization. In our 
implementation the SR-UKF is more efficient than other UKF- 
based methods, mainly because it makes use of powerful 
linear algebra techniques including the orthogonal-triangular 
decomposition and Cholesky factor updating. 

It is seen from Table m that the additional calculation 
for ‘nearPD’ is almost negligible and the computational 
complexity of UKF-GPS should also be 0(n^). For UKF- 
GPS, the number of average times that it is requires to execute 
the ‘nearPD’ algorithm in one estimation and the average 
time steps that need to execute ‘nearPD’ calculation are listed 
in Table |V] Note that in each time step ‘nearPD’ can be 
calculated before 0 or 0 in Algorithm and thus the 
number of times for executing ‘nearPD’ can be greater than 
the number of time steps involved for ‘nearPD’ calculation. 
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TABLE IV 

Time for Estimation for NPCC 48-Machine System 


Filter 

Time (second) 

three-phase fault random fault 

EKF 

42.615 

42.661 

UKF-schol 

118.580 

119.246 

UKF-k 

118.802 

119.544. 

UKF-modified 

118.806 

119.789 

UKF-AQ 

121.188 

122.274 

UKF-GPS 

119.230 

119.085 

SR-UKF 

104.733 

105.360 

TABLE V 

Average Times of Executing ‘nearPD’ and Average Number of 
Time Steps Involved in One Estimation 


Time (second) 


three-phase fault 

random fault 

Average times of 
executing ‘nearPD’ 

8.88 

8.50 

Average number of 
time steps involved 

7.34 

6.90 


VI. Conclusion 

In this paper, we introduce and compare six approaches to 
enhance the numerical stability and further the scalability of 
the unscented Kalman hlter, including the proposed UKF-GPS 
method. These methods and the extended Kalman Filter are 
tested on WSCC 3-machine system and NPCC 48-machine 
system. For WSCC system, there is no numerical stability 
problem for classic UKF, and all methods work well. However, 
for NPCC system, EKF cannot converge and UKF encounters 
numerical stability problem. Among the introduced methods, 
UKF-schol, UKF -K, and UKF-AQ can have big estimation 
errors for several estimations and UKF-AQ even diverge in 
some cases; UKF-GPS works well for almost all estimations; 
and UKF-modified and SR-UKF work very well for all esti¬ 
mations due to their better numerical stability and scalability. 

Apart from the EKF and UKF that are discussed in this 
paper, recently some other approaches have also been applied 
to dynamic state estimation, such as the extended particle 
hlter | |42l , cubature Kalman hlter | |45l , and observers | |45l , 
ph) . EKF, SR-UKF, CKF, and nonlinear observers has been 
compared for power system DSE under model uncertainty and 
malicious cyber attacks in Eg. A good comparison of EKF, 
classic UKF, ensemble Kalman hlter, and particle hlter is also 
performed in pT) . It would be valuable to more thoroughly 
compare the approaches discussed in this paper with other 
approaches in order to provide a guideline about how to choose 
the most suitable approaches for power system DSE. 
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